home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dsyev.f < prev    next >
Text File  |  1997-06-25  |  6KB  |  199 lines
  1.       SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
  2. *
  3. *  -- LAPACK driver routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          JOBZ, UPLO
  10.       INTEGER            INFO, LDA, LWORK, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  DSYEV computes all eigenvalues and, optionally, eigenvectors of a
  20. *  real symmetric matrix A.
  21. *
  22. *  Arguments
  23. *  =========
  24. *
  25. *  JOBZ    (input) CHARACTER*1
  26. *          = 'N':  Compute eigenvalues only;
  27. *          = 'V':  Compute eigenvalues and eigenvectors.
  28. *
  29. *  UPLO    (input) CHARACTER*1
  30. *          = 'U':  Upper triangle of A is stored;
  31. *          = 'L':  Lower triangle of A is stored.
  32. *
  33. *  N       (input) INTEGER
  34. *          The order of the matrix A.  N >= 0.
  35. *
  36. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
  37. *          On entry, the symmetric matrix A.  If UPLO = 'U', the
  38. *          leading N-by-N upper triangular part of A contains the
  39. *          upper triangular part of the matrix A.  If UPLO = 'L',
  40. *          the leading N-by-N lower triangular part of A contains
  41. *          the lower triangular part of the matrix A.
  42. *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
  43. *          orthonormal eigenvectors of the matrix A.
  44. *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
  45. *          or the upper triangle (if UPLO='U') of A, including the
  46. *          diagonal, is destroyed.
  47. *
  48. *  LDA     (input) INTEGER
  49. *          The leading dimension of the array A.  LDA >= max(1,N).
  50. *
  51. *  W       (output) DOUBLE PRECISION array, dimension (N)
  52. *          If INFO = 0, the eigenvalues in ascending order.
  53. *
  54. *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  55. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  56. *
  57. *  LWORK   (input) INTEGER
  58. *          The length of the array WORK.  LWORK >= max(1,3*N-1).
  59. *          For optimal efficiency, LWORK >= (NB+2)*N,
  60. *          where NB is the blocksize for DSYTRD returned by ILAENV.
  61. *
  62. *  INFO    (output) INTEGER
  63. *          = 0:  successful exit
  64. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  65. *          > 0:  if INFO = i, the algorithm failed to converge; i
  66. *                off-diagonal elements of an intermediate tridiagonal
  67. *                form did not converge to zero.
  68. *
  69. *  =====================================================================
  70. *
  71. *     .. Parameters ..
  72.       DOUBLE PRECISION   ZERO, ONE
  73.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  74. *     ..
  75. *     .. Local Scalars ..
  76.       LOGICAL            LOWER, WANTZ
  77.       INTEGER            IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
  78.      $                   LLWORK, LOPT
  79.       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
  80.      $                   SMLNUM
  81. *     ..
  82. *     .. External Functions ..
  83.       LOGICAL            LSAME
  84.       DOUBLE PRECISION   DLAMCH, DLANSY
  85.       EXTERNAL           LSAME, DLAMCH, DLANSY
  86. *     ..
  87. *     .. External Subroutines ..
  88.       EXTERNAL           DLASCL, DORGTR, DSCAL, DSTEQR, DSTERF, DSYTRD,
  89.      $                   XERBLA
  90. *     ..
  91. *     .. Intrinsic Functions ..
  92.       INTRINSIC          MAX, SQRT
  93. *     ..
  94. *     .. Executable Statements ..
  95. *
  96. *     Test the input parameters.
  97. *
  98.       WANTZ = LSAME( JOBZ, 'V' )
  99.       LOWER = LSAME( UPLO, 'L' )
  100. *
  101.       INFO = 0
  102.       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
  103.          INFO = -1
  104.       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
  105.          INFO = -2
  106.       ELSE IF( N.LT.0 ) THEN
  107.          INFO = -3
  108.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  109.          INFO = -5
  110.       ELSE IF( LWORK.LT.MAX( 1, 3*N-1 ) ) THEN
  111.          INFO = -8
  112.       END IF
  113. *
  114.       IF( INFO.NE.0 ) THEN
  115.          CALL XERBLA( 'DSYEV ', -INFO )
  116.          RETURN
  117.       END IF
  118. *
  119. *     Quick return if possible
  120. *
  121.       IF( N.EQ.0 ) THEN
  122.          WORK( 1 ) = 1
  123.          RETURN
  124.       END IF
  125. *
  126.       IF( N.EQ.1 ) THEN
  127.          W( 1 ) = A( 1, 1 )
  128.          WORK( 1 ) = 3
  129.          IF( WANTZ )
  130.      $      A( 1, 1 ) = ONE
  131.          RETURN
  132.       END IF
  133. *
  134. *     Get machine constants.
  135. *
  136.       SAFMIN = DLAMCH( 'Safe minimum' )
  137.       EPS = DLAMCH( 'Precision' )
  138.       SMLNUM = SAFMIN / EPS
  139.       BIGNUM = ONE / SMLNUM
  140.       RMIN = SQRT( SMLNUM )
  141.       RMAX = SQRT( BIGNUM )
  142. *
  143. *     Scale matrix to allowable range, if necessary.
  144. *
  145.       ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
  146.       ISCALE = 0
  147.       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
  148.          ISCALE = 1
  149.          SIGMA = RMIN / ANRM
  150.       ELSE IF( ANRM.GT.RMAX ) THEN
  151.          ISCALE = 1
  152.          SIGMA = RMAX / ANRM
  153.       END IF
  154.       IF( ISCALE.EQ.1 )
  155.      $   CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
  156. *
  157. *     Call DSYTRD to reduce symmetric matrix to tridiagonal form.
  158. *
  159.       INDE = 1
  160.       INDTAU = INDE + N
  161.       INDWRK = INDTAU + N
  162.       LLWORK = LWORK - INDWRK + 1
  163.       CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
  164.      $             WORK( INDWRK ), LLWORK, IINFO )
  165.       LOPT = 2*N + WORK( INDWRK )
  166. *
  167. *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
  168. *     DORGTR to generate the orthogonal matrix, then call DSTEQR.
  169. *
  170.       IF( .NOT.WANTZ ) THEN
  171.          CALL DSTERF( N, W, WORK( INDE ), INFO )
  172.       ELSE
  173.          CALL DORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
  174.      $                LLWORK, IINFO )
  175.          CALL DSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ),
  176.      $                INFO )
  177.       END IF
  178. *
  179. *     If matrix was scaled, then rescale eigenvalues appropriately.
  180. *
  181.       IF( ISCALE.EQ.1 ) THEN
  182.          IF( INFO.EQ.0 ) THEN
  183.             IMAX = N
  184.          ELSE
  185.             IMAX = INFO - 1
  186.          END IF
  187.          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
  188.       END IF
  189. *
  190. *     Set WORK(1) to optimal workspace size.
  191. *
  192.       WORK( 1 ) = MAX( 3*N-1, LOPT )
  193. *
  194.       RETURN
  195. *
  196. *     End of DSYEV
  197. *
  198.       END
  199.